##########################################################################################

library(data.table)
library(optparse)
library(ggplot2)
library(ggsci)

##########################################################################################
option_list <- list(
    make_option(c("--input_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    input_file <- "~/20231121_singleMuti/results/celltype_plot/peak_region/encode_enrich/region_enrich.tsv"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot/peak_region/encode_enrich"

}


###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
out_path <- opt$out_path

dir.create( out_path , recursive = T)
#dir.create( paste0(out_path , "/peakMA") , recursive = T)


###########################################################################################

dat <- unique(data.frame(fread(input_file)))

## 38基因组的长度
hg38_len <- 3088286480
options(scipen = 999)

###########################################################################################
result <- c()
for( type in unique(dat$name1) ){
    print(type)

    dat_use <- subset( dat , name1 == type )
    dat_use$P <- ""
    dat_use$ER <- ""

    for( i in 1:nrow(dat_use)){
        print(i)
        ## 关注peak在注释区域上
        a <- dat_use[i,"total_coverage1"] - dat_use[i,"uniq_coverage1"]
        ## 关注peak不在注释区域上
        b <- dat_use[i,"uniq_coverage1"]
        ## 其它区域在注释区域上
        c <- dat_use[i,"uniq_coverage2"]
        ## 其它区域不在注释区域上
        d <- hg38_len - a - b - c
        data <- matrix( c(a,b,c,d) , ncol =2 )
        chisq_result <- chisq.test(data)

        p_value <- chisq_result$p.value

        ## 富集系数
        enrichment <- (a/(a+b)) / (c/(c+d))
        dat_use$P[i] <- p_value
        dat_use$ER[i] <- enrichment
    }

    result <- rbind( result , dat_use )
}

images_name <- paste0( out_path , "/histon_enrich.tsv" )
write.table( result , images_name , row.names = F , quote = F , sep = "\t" )

###########################################################################################
## 柱状图可视化富集
result$name2 <- sapply(strsplit(result$name2 , "_merge") , "[" , 1)

col <- pal_lancet("lanonc")(9)
names(col) <- unique(result$name2)[length(unique(result$name2)):1]

for( type in unique(result$name1) ){
    print(type)

    result_plot <- subset( result , name1 == type )
    result_plot$log2_ER <- log2(as.numeric(result_plot$ER))
    result_plot$log2_ER <- ifelse(result_plot$ER < 1 , -1 , result_plot$log2_ER )

    order_histon <- result_plot[order(result_plot$log2_ER , decreasing = T),"name2"]
    result_plot$name2 <- factor( result_plot$name2 , levels = order_histon , order = T )

    p <- ggplot(result_plot , mapping = aes( name2 , log2_ER , fill = name2)) +
        geom_bar(stat='identity' , position='stack' , color = "black") +
        scale_fill_manual(values=col) + 
        ylim(-1.5 , 1.5 * max(result_plot$log2_ER)) +
        xlab("") +
        ylab("log2 fold enrichment") +
        theme(
                title =element_text(size=4, face='bold'),
                legend.position = 'none' ,
                panel.grid.major=element_blank(),
                panel.grid.minor=element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1 , size = 15  , color = 'black' , family="Helvetica" ) ,
                axis.title.y =  element_text(size = 15) ,
                axis.text.y =  element_text(size = 15 , color = 'black' , family="Helvetica") ,
                axis.ticks.length=unit(.25, "cm") ,
                axis.line = element_line(size = 0.5)
        )

    images_name <- paste0( out_path , "/histon_enrich." , type , ".pdf" )
    ggsave( images_name , p , width = 6 , height = 5 )

}